class: center, middle, inverse, title-slide # Manipulación y generación de mapas en R ## Curso: Introducción al análisis espacial y web-mapping
con Google Earth Engine y R Shiny
### MSc. José A. Lastra
Matías Olea
### Laboratorio Geo-Información y Percepción Remota ### 06/07/2021 --- background-image: url(data:image/png;base64,#logo_labgrs_color.png) background-position: center background-size:40% --- Librerías utilizadas en esta sesión ```r library(tidyverse) #manipulación ordenada library(sf) #manipulación de vectores mediante simple features library(rgdal) # manipulación de vectores base library(units) #da soporte para medición de unidades en objetos en R library(raster) #manipulación raster library(plotly) #librería para manipulación de gráficos dinámicos library(leaflet) #libreria para manipulación de mapas dinámicos ``` --- #Recordemos -- > To understand computations in R, two slogans are helpful: - Everything that exists is an object, and - Everything that happens is a function call. .footnote[John Chambers creador de R, en Advanced R (2019) 2da Edición, p. 79.] --- # Lectura de datos -- - Antes de comenzar con los análisis y con la información leeremos nuestros datos vectoriales a trabajar -- - En esta sesión usaremos las capas de nombres: *censoINE_Valparaiso2017.gpkg* y *puntosVerdesyLimpios_MMA.gpkg* ```r valpo <- read_sf('censoINE_Valparaiso2017.gpkg') pvl <- read_sf('puntosVerdesyLimpios_MMA.gpkg') ``` <img src="data:image/png;base64,#01_mapas_dinamicos_files/figure-html/unnamed-chunk-3-1.png" width="504" style="display: block; margin: auto;" /> --- # Datos INE -- <div style="border: 1px solid #ddd; padding: 5px; overflow-x: scroll; width:1000px; "><table class="table table-striped table-hover table-condensed table-responsive" style="margin-left: auto; margin-right: auto;"> <thead> <tr> <th style="text-align:right;"> OBJECTID </th> <th style="text-align:left;"> REGION </th> <th style="text-align:left;"> NOM_REGION </th> <th style="text-align:left;"> PROVINCIA </th> <th style="text-align:left;"> NOM_PROVIN </th> <th style="text-align:left;"> COMUNA </th> <th style="text-align:left;"> NOM_COMUNA </th> <th style="text-align:right;"> T_HOM </th> <th style="text-align:right;"> T_MUJ </th> <th style="text-align:right;"> T_POB </th> <th style="text-align:right;"> T_VIV </th> <th style="text-align:right;"> SUPERFICIE </th> <th style="text-align:right;"> Densidad_ </th> <th style="text-align:left;"> geom </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 10 </td> <td style="text-align:left;"> 5 </td> <td style="text-align:left;"> REGIÓN DE VALPARAÍSO </td> <td style="text-align:left;"> 51 </td> <td style="text-align:left;"> VALPARAÍSO </td> <td style="text-align:left;"> 5104 </td> <td style="text-align:left;"> JUAN FERNÁNDEZ </td> <td style="text-align:right;"> 499 </td> <td style="text-align:right;"> 427 </td> <td style="text-align:right;"> 926 </td> <td style="text-align:right;"> 445 </td> <td style="text-align:right;"> 103.9921 </td> <td style="text-align:right;"> 8.904521 </td> <td style="text-align:left;"> MULTIPOLYGON (((-595109.1 6... </td> </tr> <tr> <td style="text-align:right;"> 53 </td> <td style="text-align:left;"> 5 </td> <td style="text-align:left;"> REGIÓN DE VALPARAÍSO </td> <td style="text-align:left;"> 54 </td> <td style="text-align:left;"> PETORCA </td> <td style="text-align:left;"> 5403 </td> <td style="text-align:left;"> PAPUDO </td> <td style="text-align:right;"> 3341 </td> <td style="text-align:right;"> 3015 </td> <td style="text-align:right;"> 6356 </td> <td style="text-align:right;"> 5823 </td> <td style="text-align:right;"> 166.6183 </td> <td style="text-align:right;"> 38.147083 </td> <td style="text-align:left;"> MULTIPOLYGON (((279038.4 64... </td> </tr> <tr> <td style="text-align:right;"> 65 </td> <td style="text-align:left;"> 5 </td> <td style="text-align:left;"> REGIÓN DE VALPARAÍSO </td> <td style="text-align:left;"> 57 </td> <td style="text-align:left;"> SAN FELIPE DE ACONCAGUA </td> <td style="text-align:left;"> 5704 </td> <td style="text-align:left;"> PANQUEHUE </td> <td style="text-align:right;"> 3677 </td> <td style="text-align:right;"> 3596 </td> <td style="text-align:right;"> 7273 </td> <td style="text-align:right;"> 2514 </td> <td style="text-align:right;"> 121.0894 </td> <td style="text-align:right;"> 60.063068 </td> <td style="text-align:left;"> MULTIPOLYGON (((334441.9 63... </td> </tr> <tr> <td style="text-align:right;"> 67 </td> <td style="text-align:left;"> 5 </td> <td style="text-align:left;"> REGIÓN DE VALPARAÍSO </td> <td style="text-align:left;"> 54 </td> <td style="text-align:left;"> PETORCA </td> <td style="text-align:left;"> 5405 </td> <td style="text-align:left;"> ZAPALLAR </td> <td style="text-align:right;"> 3704 </td> <td style="text-align:right;"> 3635 </td> <td style="text-align:right;"> 7339 </td> <td style="text-align:right;"> 6961 </td> <td style="text-align:right;"> 289.0238 </td> <td style="text-align:right;"> 25.392368 </td> <td style="text-align:left;"> MULTIPOLYGON (((268575.2 64... </td> </tr> <tr> <td style="text-align:right;"> 73 </td> <td style="text-align:left;"> 5 </td> <td style="text-align:left;"> REGIÓN DE VALPARAÍSO </td> <td style="text-align:left;"> 52 </td> <td style="text-align:left;"> ISLA DE PASCUA </td> <td style="text-align:left;"> 5201 </td> <td style="text-align:left;"> ISLA DE PASCUA </td> <td style="text-align:right;"> 3819 </td> <td style="text-align:right;"> 3931 </td> <td style="text-align:right;"> 7750 </td> <td style="text-align:right;"> 3136 </td> <td style="text-align:right;"> 164.0956 </td> <td style="text-align:right;"> 47.228561 </td> <td style="text-align:left;"> MULTIPOLYGON (((-3697793 62... </td> </tr> </tbody> </table></div> --- # Datos MMA -- <div style="border: 1px solid #ddd; padding: 5px; overflow-x: scroll; width:1000px; "><table class="table table-striped table-hover table-condensed table-responsive" style="margin-left: auto; margin-right: auto;"> <thead> <tr> <th style="text-align:right;"> OBJECTID </th> <th style="text-align:left;"> type </th> <th style="text-align:left;"> address_ty </th> <th style="text-align:left;"> address_na </th> <th style="text-align:left;"> address_nu </th> <th style="text-align:right;"> region_id </th> <th style="text-align:left;"> commune_id </th> <th style="text-align:left;"> owner </th> <th style="text-align:left;"> manager </th> <th style="text-align:left;"> glass </th> <th style="text-align:left;"> paper </th> <th style="text-align:left;"> plastic </th> <th style="text-align:left;"> phone </th> <th style="text-align:left;"> geom </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 1031 </td> <td style="text-align:left;"> Puntos Verdes </td> <td style="text-align:left;"> Calle </td> <td style="text-align:left;"> Esmeralda </td> <td style="text-align:left;"> NA </td> <td style="text-align:right;"> 5 </td> <td style="text-align:left;"> Papudo </td> <td style="text-align:left;"> CRISTORO </td> <td style="text-align:left;"> CRISTORO </td> <td style="text-align:left;"> Si </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> MULTIPOINT ((270230.7 64007... </td> </tr> <tr> <td style="text-align:right;"> 951 </td> <td style="text-align:left;"> Puntos Verdes </td> <td style="text-align:left;"> Calle </td> <td style="text-align:left;"> Ema Lobos </td> <td style="text-align:left;"> NA </td> <td style="text-align:right;"> 5 </td> <td style="text-align:left;"> Panquehue </td> <td style="text-align:left;"> CRISTALERIAS CHILE </td> <td style="text-align:left;"> CRISTALERIAS CHILE </td> <td style="text-align:left;"> Si </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> MULTIPOINT ((321656.6 63689... </td> </tr> <tr> <td style="text-align:right;"> 954 </td> <td style="text-align:left;"> Puntos Verdes </td> <td style="text-align:left;"> Calle </td> <td style="text-align:left;"> Antofagasta </td> <td style="text-align:left;"> NA </td> <td style="text-align:right;"> 5 </td> <td style="text-align:left;"> Panquehue </td> <td style="text-align:left;"> CRISTALERIAS CHILE </td> <td style="text-align:left;"> CRISTALERIAS CHILE </td> <td style="text-align:left;"> Si </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> MULTIPOINT ((328910.1 63696... </td> </tr> <tr> <td style="text-align:right;"> 966 </td> <td style="text-align:left;"> Puntos Verdes </td> <td style="text-align:left;"> Pasaje </td> <td style="text-align:left;"> Francisco Bulnes </td> <td style="text-align:left;"> NA </td> <td style="text-align:right;"> 5 </td> <td style="text-align:left;"> Panquehue </td> <td style="text-align:left;"> CRISTALERIAS CHILE </td> <td style="text-align:left;"> CRISTALERIAS CHILE </td> <td style="text-align:left;"> Si </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> MULTIPOINT ((326368.8 63711... </td> </tr> <tr> <td style="text-align:right;"> 971 </td> <td style="text-align:left;"> Puntos Verdes </td> <td style="text-align:left;"> Calle </td> <td style="text-align:left;"> Maximiliano Errazuriz </td> <td style="text-align:left;"> NA </td> <td style="text-align:right;"> 5 </td> <td style="text-align:left;"> Panquehue </td> <td style="text-align:left;"> CRISTALERIAS CHILE </td> <td style="text-align:left;"> CRISTALERIAS CHILE </td> <td style="text-align:left;"> Si </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> No </td> <td style="text-align:left;"> MULTIPOINT ((326698.7 63714... </td> </tr> </tbody> </table></div> --- # Repaso manipulación básica vectores -- - Dentro de los análisis más elementales podemos distinguir el proceso de *dissolve* o simplificación de entidades -- - En este ejemplo podemos obtener el total de población total por provincia -- - *Importante*: haremos uso de nuestro *%>%* -- ```r provincias <- valpo %>% #objeto sf de entrada o tabla group_by(NOM_PROVIN) %>% # campo por el cual queremos simplificar summarize(poblacion = sum(T_POB)) # estadístico a calcular indicando el nombre en la nueva tabla ``` -- - Si no ponemos nada en la función *summarize()* solo estará en el nuevo objeto el nombre de los elementos en el campo usado en *group_by()* ```r provincias <- valpo %>% #objeto sf de entrada o tabla group_by(NOM_PROVIN) %>% # campo por el cual queremos simplificar summarize() #dejando vacío solo se agrupa sin ningún estadístico adicional ``` --- <img src="data:image/png;base64,#01_mapas_dinamicos_files/figure-html/unnamed-chunk-8-1.png" width="504" style="display: block; margin: auto;" /> --- # Cálculos en tabla -- - Además de alterar espacialmente la información, podemos realizar cálculos en la tabla de atricutos de nuestros objetos espaciales. -- - Al calcular el área, por defecto se calcula en metros cuadrados. Por lo cual, realizaremos la conversión a km~2 y fijaremos las unidades correctas con *set_units()* ```r area_km <- st_area(provincias) %>% set_units(km^2) provincias <- provincias %>% mutate(area_km = area_km) ``` <div style="border: 1px solid #ddd; padding: 5px; overflow-x: scroll; width:1000px; "><table class="table table-striped table-hover table-condensed table-responsive" style="margin-left: auto; margin-right: auto;"> <thead> <tr> <th style="text-align:left;"> NOM_PROVIN </th> <th style="text-align:right;"> poblacion </th> <th style="text-align:left;"> geom </th> <th style="text-align:right;"> area_km </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> ISLA DE PASCUA </td> <td style="text-align:right;"> 7750 </td> <td style="text-align:left;"> MULTIPOLYGON (((-3697793 62... </td> <td style="text-align:right;"> 245.9194 [km^2] </td> </tr> <tr> <td style="text-align:left;"> LOS ANDES </td> <td style="text-align:right;"> 110602 </td> <td style="text-align:left;"> POLYGON ((344401.3 6352420,... </td> <td style="text-align:right;"> 3053.1089 [km^2] </td> </tr> <tr> <td style="text-align:left;"> MARGA MARGA </td> <td style="text-align:right;"> 341893 </td> <td style="text-align:left;"> POLYGON ((311225.8 6337162,... </td> <td style="text-align:right;"> 1158.4480 [km^2] </td> </tr> <tr> <td style="text-align:left;"> PETORCA </td> <td style="text-align:right;"> 78299 </td> <td style="text-align:left;"> MULTIPOLYGON (((268328.3 64... </td> <td style="text-align:right;"> 4593.0597 [km^2] </td> </tr> <tr> <td style="text-align:left;"> QUILLOTA </td> <td style="text-align:right;"> 203277 </td> <td style="text-align:left;"> POLYGON ((282188.6 6366487,... </td> <td style="text-align:right;"> 1113.1464 [km^2] </td> </tr> <tr> <td style="text-align:left;"> SAN ANTONIO </td> <td style="text-align:right;"> 168046 </td> <td style="text-align:left;"> MULTIPOLYGON (((249147.9 63... </td> <td style="text-align:right;"> 1500.2938 [km^2] </td> </tr> <tr> <td style="text-align:left;"> SAN FELIPE DE ACONCAGUA </td> <td style="text-align:right;"> 154718 </td> <td style="text-align:left;"> POLYGON ((311885.3 6366821,... </td> <td style="text-align:right;"> 2634.5372 [km^2] </td> </tr> <tr> <td style="text-align:left;"> VALPARAÍSO </td> <td style="text-align:right;"> 751317 </td> <td style="text-align:left;"> MULTIPOLYGON (((-612210.4 7... </td> <td style="text-align:right;"> 2022.5651 [km^2] </td> </tr> </tbody> </table></div> --- # Mapa estático y dinámico -- - En sesiones anteriores vimos que podíamos emplear *plot()* o *geom_sf()* para graficar un objeto espacial ya sea raster o vector. -- - A continuación combinaremos el uso de *ggplot()* y la librería *plotly* para crear un mapa dinámico base -- - Mapa base con nuestros objetos sf y ggplot2 ```r # filtro de datos marga.marga <- valpo %>% filter(PROVINCIA =="58") ``` --- - *Importante*: Las capas se graficarán en el orden en que se pongan en el código ```r ggplot() + #canvas del gráfico geom_sf(data = marga.marga, aes(fill=T_POB))+ #sf con polígonos geom_sf(data = pvl)+ #sf con capa de puntos theme_bw() + scale_fill_viridis_c() #tema blanco y negro y paleta de colores para el fill ``` <img src="data:image/png;base64,#01_mapas_dinamicos_files/figure-html/unnamed-chunk-12-1.png" width="720" style="display: block; margin: auto;" /> -- - Recuerde que para el gráfico de objetos *sf* con ggplot siempre debe usar el argumento *data =* para evitar problemas con el código. --- # Convirtiendo a mapa dinámico -- - Un paso clave es guardar nuestro mapa en un objeto siguiendo el ejemplo -- ```r g <- ggplot() + #canvas del gráfico geom_sf(data = marga.marga, aes(fill=T_POB))+ #sf con polígonos geom_sf(data = pvl)+ #sf con capa de puntos theme_bw() + scale_fill_viridis_c() #tema blanco y negro y paleta de colores para el fill class(g) ``` ``` ## [1] "gg" "ggplot" ``` -- - Ahora convertiremos el objeto usando la función *ggplotly()* perteneciente a la librería [*plotly*](https://github.com/ropensci/plotly) ```r ggplotly(g) #convertir ggplot object a plotly y htmlwidget ``` ---
- En este mapa se puede hacer zoom, revisar valores, realizar snapshots, entre otras cosas integradas en plotly. - Podemos hacer mapas similares empleando capas rasters pero estas pueden ser muy pesadas para nuestra interfaz. - Estas limitantes nos llevan a trabajar con [**leaflet**](https://github.com/rstudio/leaflet) --- # Mapas dinámicos con Leaflet -- - Leaflet es una de las librerías abiertas más populares para trabajar con mapas interactivos. -- - En R disponemos de la librería *leaflet* para poder realizar los mapas dinámicos, empleando vectores (shapefiles, geopackage, geojson), capas raster (TIF, GRD, BIL, etc) y otros servicios. -- - Realicemos el código para nuestro mapa base ```r leaflet() %>% addTiles() #crea el canvas y añade tiles base de openstreet map ``` ---
- Si queremos cambiar el mapa base, podemos emplear la función *addProviderTiles()* y emplear entre caracteres algún proveedor disponible [[Ver más](http://leaflet-extras.github.io/leaflet-providers/preview/index.html)] - Una lista de los proveedores puede ser obtenida llamando al objeto *providers* --- # Añadiendo nuestras capas -- - Primero debemos reproyectar nuestra información para que opere en leaflet correctamente -- - leaflet emplea el código EPSG 4326, por lo que debemos asegurar integridad en la información a cargar. -- - Acá utilizaremos la función *st_cast()* para pasar nuestro archivo a tipo puntos ya que, hasta la fecha no hay soporte para tipo multipuntos en leaflet. -- - *Importante*: Siempre considere revisar bien las características de su capa de puntos y verifique que no haya repitición o inconsistencias en su capa. Incluso si viene desde fuentes oficiales. -- ```r valpo_latlong <- valpo %>% st_transform(4326) pvl_latlong <- pvl %>% st_cast('POINT') %>% st_transform(4326) leaflet() %>% addTiles() %>% #canvas del mapa addPolygons(data = valpo_latlong) %>% #agregando polígonos addCircleMarkers(data = pvl_latlong,opacity = 1) #agreganfo puntos ``` -- - El orden de las capas es importante, pero ese orden se puede cambiar al habilitar la activación o desactivación de capas de mapa. -- - Podemos agregar y controlar variados elementos de nuestras capas espaciales -- - Cada uno lo iremos viendo a medida que avancemos en el curso ---
--- # Leyendas -- - El control de la leyenda y la visualización es importante para comunicar la información. -- - A continuación revisaremos como configurar la leyenda de nuestros datos de población. -- - Para buscar paletas de colores crearemos en este caso un vector de colores HTML usando la función *c()* -- - Usaremos *stroke* para controlar la presencia de la línea y *weight* el ancho, *fillColor* para aplicar nuestra paleta y *markerClusterOptions()* para adaptar la visualización de los puntos verdes según la escala -- ```r # paleta de colores pal <- colorBin(palette = c('#ffffb2','#fed976','#feb24c', '#fd8d3c','#fc4e2a','#e31a1c','#b10026'), domain = valpo_latlong$T_POB, bins = 7) #Mapa leaflet() %>% addTiles() %>% addPolygons(data = valpo_latlong,stroke = T,weight = 0.3,color='black', fillOpacity = 0.8,fillColor = ~pal(T_POB)) %>% addCircleMarkers(data = pvl_latlong,opacity = 1,clusterOptions = markerClusterOptions()) %>% addLegend(pal = pal,values = valpo_latlong$T_POB,title = 'Población Total Censo 2017') ``` -- - Hay diferentes funciones para crear las paletas *colorNumeric*, *colorQuantile*, *colorFactor*. -- - Puede además, disponer de paletas pre-armadas con paquetes como *viridis* o *RColorBrewer* ---
--- #Controlando las capas en el mapa -- - Para que aparezca la información de nuestras capas y podamos activarlas o desactivarlas, usaremos el argumento *group* para asignarles un nombre distintivo y la función *addLayersControl()* para que aparezca el cuadro de control de capas. -- - Además cambiaremos la posición de la leyenda y los controles de capa -- ```r leaflet() %>% addTiles() %>% addPolygons(data = valpo_latlong,stroke = T,weight = 0.3,color='black', fillOpacity = 0.8,fillColor = ~pal(T_POB), group = 'Censo INE, Comunas') %>% #asignamos nombre addCircleMarkers(data = pvl_latlong,opacity = 1, clusterOptions = markerClusterOptions(),group = 'Puntos Verdes, MMA') %>% addLegend(pal = pal,values = valpo_latlong$T_POB, title = 'Población Total Censo 2017',position ="bottomleft") %>% addLayersControl(overlayGroups = c('Censo INE, Comunas','Puntos Verdes, MMA'), #nombres deben ser iguales position ="bottomleft") ``` -- - *Importante*: Varias capas pueden estar controladas por el mismo nombre de grupo y la legenda puede también disponer de un nombre de grupo para activarla y desactivarla. -- - Además podemos añadir varios proveedores diferentes y permitir cambiar entre los mapas base ---
--- - Ahora, terminemos agregando más proveedores usando también el argumento *baseGroups* -- ```r leaflet() %>% addProviderTiles(providers$OpenStreetMap,group = 'OSM') %>% #proveedor 1 addProviderTiles(providers$Stamen,group = 'Stamen') %>% #proveedor 2 addProviderTiles(providers$Esri.WorldImagery,group = 'Esri') %>% #proveedor 3 addPolygons(data = valpo_latlong,stroke = T,weight = 0.3,color='black', fillOpacity = 0.8,fillColor = ~pal(T_POB), group = 'Censo INE, Comunas') %>% addCircleMarkers(data = pvl_latlong,opacity = 1, clusterOptions = markerClusterOptions(),group = 'Puntos Verdes, MMA') %>% addLegend(pal = pal,values = valpo_latlong$T_POB, title = 'Población Total Censo 2017',position ="bottomleft") %>% addLayersControl(overlayGroups = c('Censo INE, Comunas','Puntos Verdes, MMA'), baseGroups = c('OSM','Stamen','Esri'),#control de proveedores position ="bottomleft") ``` ---
--- background-image: url(data:image/png;base64,#logo_labgrs_color.png) background-position: center background-size:40%